Anomalous Heat Conduction in Three-Dimensional Nonlinear Lattices 



Hayato ShibjQ 

Department of Physics, Kyoto University, Kyoto 606-8502, Japan 



00 

o 

O 

(N 



(N 



c3 



I 

C 

o 
o 



> 
(N 

(N 

o 



X 



Nobuyasu ItcQ 

Department of Applied Physics, University of Tokyo, Tokyo 113-8656, 

(Dated: May 12, 2008) 



Japan 



Heat conduction in three-dimenisional nonlinear lattice models is studied using nonequilibrium 
molecular dynamics simulations. We employ the Fermi-Pasta-Ulam-/3 model, in which nonlinearity 
exists in the interaction of the biquadratic form. It is confirmed that the thermal conductivity, 
the ratio of the energy flux to the temperature gradient, diverges with increasing system size up 
to 128 x 128 x 256 lattice sites. This size corresponds to nanoscopic to mesoscopic scales of ap- 
proximately lOOnm. From these results, we conjecture that the energy transport in insulators with 
perfect crystalline order exhibits anomalous behavior. The effects of the lattice structure, random 
impurities, and the natural length in interactions are also examined. We find that fee lattices display 
stronger divergence than simple cubic lattices. When impurity sites of infinitely large mass, which 
are thus fixed, are randomly distributed, such divergence vanishes. 



I. INTRODUCTION 

It is a long-standing problem to reproduce irreversible 
heat conduction phenomena, which are described by the 
Fourier law, J — — kVT, where k denotes the heat con- 
ductivity, on the basis of time-reversible microscopic dy- 
namics. Recent studies have succeeded in reproducing 
the Fourier law [J, H, i, 0], and other linear nonequilib- 
rium transport phenomenally ■ However, there remains 
an unsolved problem. The studies cited above mainly 
used particle systems, such as hard spheres and Lennard- 
Jones particles, for example, but historically, nonlinear 
lattice systems have also been studied@, H, H, E3, EH, 03 • 
Shimada et al. showed that the Fourier law is reproduced 
in a three-dimensional (3D) polymer-like lattice with 
an Fermi-Pasta-Ulam-/? (FPU-/?) type interaction but 
Shiba et al. observed the divergence of thermal conduc- 
tivity for 3 and 4D FPU-/? lattices [Hj]. The purpose of 
this article is to study the energy transport behavior in 
nonlinear lattice systems. 

Let us start with a brief review of the study of heat 
conduction in nonlinear lattices. It is well-known that 
no temperature gradients are formed if we use integrable 
systems, such as harmonic chains [l4j. To realize linear 
temperature gradients, we need features that give rise to 
thermalization processes. One possibility for generating 
such a feature is represented by nonintegrable nonlinear 
interactions ,15]. Nonintegrability results in macroscop- 
ically irreversible processes, even in conserved systems, 
even though the basic equation is reversible. 

In recent years, heat conduction in ID nonintegrable 
chains has been widely investigated^- In these studies, 
it was found that the heat conductivity of nonlinear lat- 
tices depends on the system size as k ~ N a , where N 



denotes the length of the system, when the total momen- 
tum of the system is conserved. Although the value of 
the exponent a seems to depend on the model, it is ap- 
proximately 0.4 for several models, for example, FPU-/? 
lattices^, FPU-a lattices diatomic Toda lattices (lfj|. 
and ID binary hard-sphere gases [IB, G3] exhibit similar 
power-law behavior. This anomalous behavior originates 
in the power-law decay of the energy-flux autocorrelation 
function, C(t) = (J(t) ■ J(0)), which is the integrand of 
the Green-Kubo formula[3 L3J, HI , 
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where J it) is the total heat flux and V is the volume 
of the system. In ID nonlinear lattices, V = N is the 
chain length. The autocorrelation functions of systems 
with a ~ 0.4 decay approximately as C(t) ~ t~° with 
/? ~ 0.6. Phenomenological explanations for this de- 
cay exponent have been proposed using mode-coupling 
theory [2l|, the renormalization group analysis of ID fluc- 
tuating hydrodynamics [H, HI], and kinetic theory [24}. 

The same types of divergence of k have been ob- 
served in ID and 2D fluid systems. For simple classi- 
cal fluids, it has been well established that the autocor- 
relation function exhibits power-law decay called "long- 
time tail" [H, HE H3, HI, Hi, El, which takes the form 
C(t) - t~ d / 2 for d> 2. This was discovered by Alder 
and Wainwright [3ll. |32| using computer simulations, and 
the study of long-time tail is still ongoing [33l |34|. In 3D 
systems, the conductivity converges in the limit N — ► oo, 
as demonstrated by the Green-Kubo formula, eq. (JTJ) . 
In 2D systems, it diverges as ln/V. Such dimensionality 
dependence has been confirmed by computer simulations 
with hard-particle systems Lennard- Jones fluids [|, 
and other systems. The system dimensionality plays a 
key role in the behavior of transport coefficents in fluid 
systems. 

How docs the conductivity of nonlinear lattice systems 
depend on their dimensionality? One may naively expect 
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a dependence that is similar to that in fluids, assuming 
that a phonon gas can be described by hydro-dynamic 
equations of motion. But this appears to be a too naive 
assumption. Recently it was shown that the conductivity 
of a 3D extension of FPU-/? lattices displays divergence 
and that this divergence is consistent with the power-law 
decay of the autocorrelation function [13j. 

Although the derivation of the long-time tail from 
the linearized hydrodynamic eq uations using the mode- 
coupling hypothesis [12, [H, H3 [H, [29|, [3(| implies its 
validity for mesoscopic to macroscopic scales, but not 
for microscopic scales, particle systems also possess long- 
time tails on microscopic scales. In simple particle sys- 
tems, distances on the order of ten times the mean- 
free path are sufficient to observe the characteristic size 
dependence of the conductivity [J, [1, In nonlinear 
lattices[13], however, it has been found that anoma- 
lous divergence continues up to systems of at least 
64 x 64 x 128 lattice sites. This size corresponds to meso- 
scopic systems. This means that the energy transport 
in nanoscopic to mesoscopic crystals is more complicated 
than that in fluid systems. From these findings, it is 
clear that further study of this problem is interesting 
both theoretically and with regard to nano technology. 
Such anomalous features should contribute to explana- 
tion of the anomalous features in heat transport observed 
in experiments and numerical simulations of nano scale 
systems such as carbon nanotubes[35l l36l |37| . 

Our model and simulation method are presented in the 
next section. The robust divergence of the thermal con- 
ductivity in 3D nonlinear lattices is demonstrated in the 
following two sections: In §3, a simple cubic lattice is 
treated, and in §4, an fee lattice is treated. Such diver- 
gence is also shown to exist in a system with a nonlinear 
interaction possessing a natural length scale in §5. In §6, 
the effect of the randomness of the mass is studied. In 
§7, the results for a 2D nonlinear lattice are given. The 
last section contains a summary and conclusion. 



II. MODEL AND SIMULATION METHOD 

The first model we study is a simple 3D extension of 
FPU-/3 lattices [13j. The model Hamiltonian is 
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where the 3D vectors pi and r% denote the momentum 
and displacement of a particle at lattice point i, respec- 
tively, and the mass m is taken to be unity. The sum- 
mation over is a sum over nearest neighbors, k and 
g are parameters indicating the strength of interactions 
between these nearest-neighbor particles, and they are 
fixed as k = 1.0 and g = 0.1 in this paper. 

We point out here that the spatial degrees of freedom 
of these dynamical models consist of the displacement 
vector Ti = qi — q®, where qi represents the real spatial 



coordinates of the particles, and q° represents the equi- 
librium positions of the particles, which are at simple 
cubic lattice points. In other words, this model system 
possesses a crystal structure. In this model the longitu- 
dinal and transverse modes are treated identically, and 
thus have the same dispersion relation in the harmonic 
limit. 

The system size is denoted by N x x N y x N z . We ex- 
press the system size in terms of the number of particles 
(N x ,N y , N z ), not the lengths (L x , L y , L z ). Thus, the sys- 
tem size is dimensionless N c = L c /a (c = x,y,z). Here, 
a is the lattice spacing constant, which does not appear 
in the Hamiltonian given by eq. ([2]). The only character- 
istic length scale that appears in eq. @ is \fk~fg, which 
is a typical length of phonon-phonon interactions. 

We use periodic boundary conditions in the x- and y- 
directions. The particles on both ends in the z-direction 
are attached to rigid walls, which are separated by one 
lattice space from the walls. These particles interact with 
the wall through a potential that is identical to that for 
particle-particle interactions, and their local temperature 
is controlled by the Nose-Hoover method [3^|. Thus, the 
equations of motion are modified for these particles as 
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Here, d denotes the Nose-Hoover thermostat variables, 
which obey 
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where T denotes the temperatures of the heat baths, 
which are set to Tl on the left and Tr on the right. 
(Here, i z increases to the right.) The difference between 
Tl and Tr, drives the energy transport in the system. 
The quantity Q is the relaxation time of the heat baths, 
which is set to unity. 

For particles in the bulk, the equations of motion are 
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m ' 
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Particle dynamics simulations were used to study this 
system. The initial displacements were zero for all par- 
ticles (i.e., Ti = 0) and the initial momenta pi were 
randomly assigned for each particle. Starting from this 
initial state, initial relaxation steps were discarded. Af- 
ter the system reached a steady state, the temperature 
distribution, energy flux, and thermal conductivity were 
computed. 

a. Temperature Using the Virial theorem, we define 
the local temperature T(i) = T(i x ,i y ,i z ) of a particle 
on lattice point i as the long-time average of the kinetic 
energy. (Here, i x , i y . and i z are the labels of the lattice 
sites, with 1 < i x < N x 
N z .) 



1 < i y < N y , and 1 < i z < 
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To obtain better accuracy, we averaged the temperature 
over the N x x N y particles in the same cross section, and 
thus define T(i z ) as 



-k B T(i z 



N X Ny 



b. Heat flux In this paper, we refer to energy flux 
as the heat flux, because the system itself is a thermody- 
namic system that relaxes to equilibrium due to its non- 
integrability. Because the heat flux is conserved in the 
bulk region, there are two ways to measure it. One way 
is to sum up the local contributions to the energy flow 
due to the interactions. The microscopic energy transfer 
from site i to site k is given by 
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where Vik is the interaction potential between sites i and 
k. By performing this summation, we could obtain the 
total heat flow density as (j z ) = N^ 1 Y^u &\ ji^k- How- 
ever, in this study, we used an alternative method; we 
calculated the work done by the heat baths as 
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Here, the summation is over all the particles to which 
the left (or right) heat baths are attached. We have con- 
firmed that these two methods yield the same value of 
the heat flow to at least two significant figures. 

c. Thermal conductivity In the simulations carried 
out for this study, the temperature distribution was 
found to be linear in the system and no temperature gap 
was observed. Therefore, the thermal conductivity k(N z ) 
can be estimated by 



(jz)N z 
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For the numerical integration of the equations of mo- 
tion, we used a Stormer-Verlet difference scheme, given 
by 
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Here, n represents the number of time steps. This 
scheme exhibits second-order convergence with At. For 
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FIG. 1: Temperature profile for the FPU-/3 model with a 
simple cubic lattice. The sequences represent the results for 
lattices with N z =48, 96, and 192, from top to bottom. The 
horizontal axis represents the position along the z-direction, 
rescaled by the system size N z . The local temperature is 
averaged over a cross-sectional cut in the a;j/-plane. The 3<r 
width is smaller than the symbols. 



the bulk particles, without temperature control, this 
scheme reduces to the leap-frog integrator, which is sym- 
plectic. Therefore, the integration of the bulk area is 
numerically stable for all At satisfying At < 0.02. How- 
ever, this scheme is not symplectic for the thermostatted 
particles. The time step At is taken to be small enough 
for the integration of these particles connected to the heat 
bath to be stable. Throughout this paper, At is set to 
0.02. 



III. SIMPLE CUBIC LATTICE 

The results for the simple cubic FPU-/? lattice are given 
in this section up to N z — 256. (A previous study gave 
results up to 128[13|.) For each data point, the results are 
averaged over simulations starting from five independent 
initial conditions. 

First, we present the results for the temperatures 
(T L ,T R ) = (20.0,10.0). For this case, the spatial tem- 
perature profile is plotted in Fig. [TJ for systems with 
N z — 48, 96, and 192 with the positions of the parti- 
cles rescaled by the system size N z . No boundary gap 
is observed in this temperature profile, and the results 
T(i z /N z ) for all values of N z nearly coincide. We thus 
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FIG. 2: (a) System size dependence of the thermal conductivity for the simple cubic FPU-/3 model plotted on a semi-log scale, 
(b) The same data plotted on a log-log scale. A power-law fitting was carried out for the data with N z > 96. The result is 
K ~ N° - 221 ( 4 ) ) which is plotted by on a broken line. 



conclude that we can estimate the thermal conductivity 
using eq. (fTO]). The system size dependence of the esti- 
mated values of k(N z ) are plotted in Fig. [2] on semi-log 
and log-log scales, up to the system size 128 x 128 x 256. 
These results for k(N z ) suggest logarithmic divergence, 
although power-law divergence is not excluded. When a 
power law is fitted to the data for 96 < N z < 256, we 
obtain k ~ n®- 221 ^\ Such a divergence was previously 
observed up to N z = 128[l3|], and it is now confirmed 
up to N z = 256. Note that the results show no sign of 
convergence. 



The behavior of k(N z ) when the heat bath tempera- 
tures are changed is shown in Fig. [3] For the same N z , 
thermal conductivity becomes lower for higher (T/^Tr). 
This is expected because thermalization is enhanced 
when the temperature is high because of the stronger 
nonlinear interactions. However, we observe that the 
divergence of k(N z ) with increasing system size is ob- 
served even at higher temperatures such as (Tl,Tr) = 
(40.0,20.0). 



After beginning the simulation with the initial con- 
ditions described in §2, we waited for the time period 
t w ~ 10 5 . By this time, the system had relaxed to a 
nonequilibrium steady state, where the amount of energy 
flow per time was stationary. For the 128 x 128 x 256 sys- 
tem, which is the largest system considered in the present 
study, and includes about 4 x 10 6 particles, the total 
number of simulation steps multiplied by the number of 
particles is of the order of 10 13 . For the 128 x 128 x 256 
system, approximately 20h of CPU time was required to 
numerically calculate one sample, using an SX8 cluster 
with a peak performance of 128 GFlops. 



IV. FCC LATTICE 

We also studied the thermal conductivity of an fee sys- 
tem in order to determine its dependence on the lattice 
structure. The Hamiltonian of this model is also that 
given in eq. but in this case, each particle inter- 

acts with its twelve nearest neighbors. Lattices of sizes 
N x N x (N + 1) are used, and each lattice point is again 
denoted by (i x , i y ,i z ). In these simulations, particles only 
occupy sites where (i x + i y + i z ) is even. Therefore, the 
total number of particles is N 2 (N + l)/2. Nose-Hoover 
heat baths were attached to the 1st and (N + l)th layers 
in the z-direction. In these layers, the particles are con- 
nected to the walls with the same interactions as those in 
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FIG. 3: System size dependence of the thermal conductivity 
for the temperatures (T L ,T R ) = (10.0, 5.0), (20.0, 10.0), and 
(40.0, 20.0). The data for (T L ,T R ) = (20.0, 10.0) are the same 
as those in Fig. [2] 
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FIG. 4: Size dependence of thermal conductivity of an fee FPU-/3 nonlinear lattice on (a) semi-log and (b) log-log scales. The 
fit of the result for N > 64 is plotted with the dotted line, which represents k ~ jV 0,51(1) . 



the system discussed above (fixed boundary condition). 
Periodic boundary conditions were adopted in the x- and 
y-directions. Similarly to the simple cubic case, It takes 
a time period of less than 10 5 for the system to reach a 
nonequilibrium steady state. 

The estimated values of the thermal conductivity are 
shown in Fig. Q]up to the system size of 192 x 192 x 193. 
The temperatures of the heat baths here are (Tl,Tr) = 
(20.0, 10.0). In this figure, power-law divergence is clearly 
observed, and in this case, the possibility of logarithmic 
divergence is excluded. Fitting in the region 64 < N < 
192 yields k(N) ~ iV ' 51 ^). This value of the exponent 
is larger than that for the simple cubic lattice 0.221(4). 
This suggests that the power-law exponent of the diver- 
gence is not determined only by the lattice dimension- 
ality. It is interesting that this divergence seems to be 
stronger than that in ID models, for which the exponent 
was found to be approximately 0.40. 



V. NATURAL LENGTH 

Using the Hamiltonian in eq. @, we cannot realize 
a system in which longitudinal and transverse modes are 
treated differently. However, in real solids, the dispersion 
relations for these modes generally differ, and they can 
have very different structures. To include such effects, 
we need to consider the interactions among the particles 
more carefully. 

Here, we consider the inclusion of the natural length 
of the interaction as a first step, and we consider the 
modified Hamiltonian 
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where dq^ = \qi — Qj \ — lo: an d 'o 1S the natural length of 
the interaction. A system comprising a simple cubic lat- 



tice was investigated using this Hamiltonian by Shimada 
et ai, and they observed normal heat conduction in three 
dimensions [l|. However, this model has the flaw that its 
transverse modes are softened, and as a result, the crystal 
structure was not maintained in their simulations. 

In order to realize a system that maintains a crystal 
structure in thermal equilibrium employing such a class 
of Hamiltonians, we need to employ a stable structure. 
One possibility is the fee lattice, which can maintain its 
structure with only nearest-neighbor interactions. In this 
section, we present simulation results for systems com- 
prising fee lattices, employing the Hamiltonian in eq. 
(fTTIl . Here, we set the parameters as 



k = 1.0, g = 0.1, l Q = 50.0 x V2. 



(12) 




FIG. 5: Temperature profiles of the fee nonlinear lattice whose 
Hamilitonian is given by eq. (|11[) . Results for the sizes TV = 64 
and 128 are displayed. Their curves are indistinguishable in 
this plot. 
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FIG. 6: System size dependence of the thermal conductivity for an fee lattice whose Hamiltonian is given by eq. (Ilip . On the 
left is a plot on a semi-log scale, and on the right is a plot of the same data on a log-log scale. In the right figure, the dotted 
line represents the result of a power-law fitting for the region N > 48 giving the result k ~ jV ' 295 ' 4 '. 



The lattice structure here is the same as that in the pre- 
vious sections. Thus, the heat flux flows along one of 
the face-centered directions, and the lattice constant is 
V2lo — 100.0. The temperatures of the heat baths were 
also set to (T L ,T R ) = (20.0, 10.0) in this case. 

The temperature profiles are plotted in Fig. [5] We 
again observe that the temperature gaps near the walls 
are small and that the results T(i z /N z ) are the same 
for different system sizes N. We can therefore obtain the 
system size dependence of thermal conductivity using eq. 
([10]). As in the previous case, we consider Nx N x (N+ 1) 
lattices with N 2 (N + l)/2 particles. 

The system size dependence of thermal conductivity 
is plotted in Fig. [5] on both semi-log and log-log scales. 
It is clearly seen that the divergence is not logarithmic 
but of a power-law type. Fitting the data for the region 
N > 48, we obtain a tentative value for the power-law 
exponent of a ~ 0.295(4). 

This result suggests that the anomalous behavior we 
have observed is not peculiar to a Hamiltonian of the 
form of eq. ([2]), and we conjecture that it is a general 
property of insulating solids. So far, all the results we 
have obtained for nonlinear lattice models with perfect 
crystalline order are consistent with anomalous thermal 
conductivity. It is a future problem to investigate the 
extent to which this anomalous behavior can be observed. 



VI. DISORDER 

A. Random mass 

The purpose of this paper is to show that the diver- 
gence of thermal conductivity is a robust property of 
nonlinear lattices, even in 3D cases. It has been a long- 
standing belief that the following play important roles 
in thermal and transport phenomena: (i) nonlinearity in 



the interactions, (ii) the dimension of the system, and 
(hi) impurities or disorders. In the previous sections, we 
showed that anomalous behavior exists if we include only 
(i) and (ii). In this section, we investigate some systems 
with disorder. 

In ID cases, disordered harmonic chains have been 
investigated for many years [1^, 0, For such 

systems driven by a Langevin thermostat, it has been 
proved that there is a unique steady state with k(N) ~ 

Li et al. investigated heat transport in disordered FPU 
chains using Nose- Hoover thermostats 43] . They investi- 
gated heat transport at various temperatures and found 
that for sufficiently high temperatures, the thermal con- 
ductivity diverges as k(N) ~ TV 43 . Although in the 
case of Nose-Hoover thermostats, at low temperatures a 
unique noncquilibrium steady state might not exist, when 
the temperature T is much higher than 0.5, there should 
be no such problem. To eliminate this possible problem, 
therefore, we only investigated 3D systems at high tem- 
peratures. The model we investigated has the following 
Hamiltonian: 
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The only difference between this Hamiltonian and that 
appearing in eq. ^ is that, in the above form, the masses 
of the particles m t vary among lattice sites. This allows 
us to include disorder in the form of random masses. We 
chose each m; as a random number given by 



to + X(Ri - 0.5), 



(14) 



where A is a parameter adjusting the amplitude of ran- 
domness, and Ri is a random number distributed uni- 
formly on the interval of [0, 1). We set the average mass 
Too to unity. 
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FIG. 7: (a) System size dependence of thermal conductivity of 3D disordered FPU-/3 lattices with the Hamiltonian given by 
eq. (|13[) . plotted on semi-log scale. The cases of two strengths of the disorder parameter A = 0.2 and 0.4 are investigated, (b) 
Plot of the same data on log-log scale. 



Using the Hamiltonian in eq. (|13p , we carried out sim- 
ulations for a system with a simple cubic lattice. The 
setup for the simulation was similar to that for the sim- 
ulations discussed in §3, with N x : N y : N z = 1 : 1 : 2. 
The averages of the thermal conductivities obtained in 
5 mass configurations were computed for various system 
sizes and two values of A. The system size dependence 
of the thermal conductivity is plotted in Fig. [71 It is ob- 
served that the thermal conductivity is lower when the 
effect of the disorder is stronger and that, here too, ther- 
mal conductivity diverges with increasing system size. 
We conclude that disorder does not destroy the anoma- 
lous thermal conductivity. 



B. Randomly fixed sites 

Because we could find no indication of normal heat 
conduction in the simulations considered so far, we also 
studied one extreme situation when the disorder effect 
should be very strong, that in which the system possesses 
impurities of infinite mass. Although the form of the 
Hamiltonian in this case is also given by eq. (113|) . the 
mass was chosen as 

nii = 1 or oo, (15) 

with a certain fraction of the masses set to infinity. In 
other words, we fixed a certain fraction of the particles 
at their mechanical equilibrium position with = 0. We 
selected the fixed particles randomly from all the parti- 
cles that are not connected to the heat baths with pro- 
portions of 10% and 20%. These values are both below 
the critical percentage above which the paths of the heat 
flow are completely blocked. 



The system size dependence of the thermal conductiv- 
ity is plotted in Fig. [5J It is seen that although the model 
with 10% fixed sites exhibits some divergence, the ther- 
mal conductivity of the model with 20% fixed sites clearly 
converges to a finite value as N z is increased. These re- 
sults suggest that defects in crystalline solids play an 
important role in thermal conductivity behavior. More- 
over, the conductivity clearly depends on the density of 
the fixed sites. 

We now remark on past works investigating the 2D har- 
monic model with missing bond defect s[44j (not missing 
particles), in which each particle has a scalar (not vec- 
tor) dynamical variable, (pi, qi). This model exhibits con- 
vergent thermal conductivity when there are sufficiently 
many missing bond defects, similarly to the system stud- 
ied here. More accurate, quantitative, and detailed in- 
vestigations of the role of the various type of defects will 
be studied in the future. 



VII. TWO-DIMENSIONAL LATTICE 

Here, we comment on some investigations of 2D sys- 
tems. We have already presented evidence that a 3D 
nonlinear lattice system generally exhibits power-law de- 
pendence on its system size. Because it is unlikely that 
the thermal conductivity of a 2D system has weaker di- 
vergence than that of a three-dimensional system, we 
conjecture that 2D nonlinear lattice systems also exhibit 
power-law divergence of the thermal conductivity. 

There have not been many studies on 2D systems. Al- 
though it is not yet proven, it is widely believed that the 
thermal conductivity for such systems diverges logarith- 
mically as the size increases 0, l45l l46j. 

We simulated a system consisting of a simple square 2D 
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FIG. 8: (a) System size dependence of the thermal conductivity for a model with fixed sites. The fraction of fixed sites is 10%. 
(b) As (a), but with 20% fixed sites. It is seen that the conductivity approaches a constant value with increasing N z . 




FIG. 9: Temperature profiles for 2D FPU-/3 lattices with 
N x : N y = 1:2. The sequences represent the results for 
the sizes N y = 192, 384, and 768. The temperatures of the 
heat baths at both ends were fixed to {T L ,T R ) = (20.0, 10.0). 
The horizontal axis represents the position in the y-direction, 
scaled by the system size N y , and the vertical axis represents 
the local temperature. 



FPU-/3 lattice whose Hamiltonian is of the same form as 
eq. ([2]). We fixed the aspect ratio to N x : N y = 1 : 2, 
where the y-axis is taken to be the direction of heat flow. 
For the size 384 x 768, we waited for approximately t w ~ 
5 x 10 5 for the system to reach the nonequilibrium steady 
state. 

The temperature profile of the steady state is plotted in 
Fig. [5J It is seen that there are no temperature gaps near 
the walls. Thus, we can define the thermal conductivity 
by eq. ([TO]). 

The system size dependence of the thermal conductiv- 
ity is plotted in Fig. [IJJ There, we see that the conduc- 
tivity exhibits power-law divergence. Because the values 
of N y used here are much larger than the values of N z 
used in the three-dimensional systems that we studied, 



FIG. 10: System size dependence of the thermal conductivity 
for 2D FPU-/3 lattices plotted on a log-log scale. The dashed 
line represents the result of a power-law fitting in the region 
N y > 128, yielding the result k(N v ) ~ Ny 267(5) . 



we can clearly distinguish the behavior here from log- 
arithmic divergence. The data presented here are suf- 
ficient to conclude that the divergence is stronger than 
conventional logarithmic divergence. Our estimation of 
the power-law exponent is a ~ 0.268(3), obtained by fit- 
ting in the region N y > 128. 



VIII. SUMMARY AND CONCLUSION 

Heat conduction in FPU-/3 lattices was studied using 
nonequilibrium molecular dynamics simulations. The di- 
vergence of the thermal conductivity was observed in sim- 
ple cubic lattices up to the size 128 x 128 x 256. Such 
divergence was also observed in fee lattice systems up 
to 192 x 192 x 193, in fee lattice systems with a natu- 
ral length included in the nonlinear interaction, and in 
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simple cubic lattice systems with randomly distributed 
particle masses with a variance of 20% (A = 0.4). Simi- 
lar divergence for square lattice systems was also found 
for sizes up to 384 x 768. These divergences are character- 
ized by power-law dependences on the system size, and 
the power exponents do not appear to be unique, rang- 
ing from (logarithmic) to 0.22 for simple cubic lattice 
systems, and being approximately 0.5 for fee lattice sys- 
tems. Convergence was observed in simple cubic lattice 
systems with randomly fixed particles. 

The present results imply that the long-time-tail decay 
exponent — d/2 of the energy flux autocorrelation func- 
tion is not observed in the present systems with diverg- 
ing thermal conductivity, although it has been observed 
in particle systems with an order-of-magnitude smaller 
number of degrees of freedom[l|,|2j,|3j]. Of course, the size 
limitation is always an issue to consider in finite-system 
analysis. Decay with an exponent of —3/2 may appear 
in larger 3D nonlinear lattices, and thermal conductivity 
may converge. For this reason, it would be worthwhile 
carrying out simulations on larger lattices using more 
powerful computers. However, it should be remarked 
here that the 128 x 128 x 256 lattice is already mesoscopic 
to macroscopic in the present technological sense: this 



lattice corresponds to insulators of size 50 x 50 x lOOnm 3 
when the lattice constant is assumed to be of the scale 
of diamonds, for example. Experimental observations of 
heat transport in crystals of this scale should be carried 
out. 

There remain several open problems in addition to the 
system size problem. One is to elucidate the direction 
dependence of the heat transport; in the present study, 
the heat flow was always in the [001] direction. Another is 
to make the interaction potential function more realistic. 
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